Novel finite-differencing techniques for numerical relativity: 
application to black hole excision 
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We use rigorous techniques from numerical analysis of hyperbolic equations in bounded domains 
to construct stable finite-difference schemes for Numerical Relativity, in particular for their use in 
black hole excision. As an application, we present 3D simulations of a scalar field propagating in a 
Schwarzschild black hole background. 
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The numerical implementation of Einstein's equations 
represents a daunting task. The involved nature of the 
equations themselves, and a number of technical diffi- 
culties (related to the necessarily finite computational 
domain, the limitations in relative computational power 
and the need to deal with singularities, constraints and 
gauge freedom) imply a significant challenge. 

Numerical solutions of Einstein's equations involve 
solving a nonlinear set of partial differential equations 
on a bounded domain, and formally constitute an initial- 
boundary- value problem (IBVP). Constructing stable 
and long term well behaved numerical approximations for 
such systems with boundaries is highly non-trivial. Here 
the term numerical stability is used in the sense that is 
equivalent, through Lax's theorem, to convergence (that 
is, the property that the numerical solution will approach 
the continuum one when resolution is increased). Deli- 
cate complications arise due to corner and edges at outer 
(and possible inner) boundaries, all of which introduce 
subtleties for a stable implementation. 

Recently, however, several sophisticated tools of rigor- 
ous numerical analysis have been developed for systemat- 
ically constructing stable numerical schemes for IBVP's. 
They employ a discrete form of well-posedness and thus 
are stable by construction, at least for linear systems. At 
this time, these tools have essentially not been used by 
the numerical relativity community. The purpose of this 
paper is to outline their use, in particular for black hole 
excision. 

An IBVP consists of three ingredients: a partial dif- 
ferential equation, initial and boundary data. It is well- 
posed if a solution exists, is unique, and depends con- 
tinuously on the initial and boundary data. It is well 
known that stable finite difference schemes approximat- 
ing an IBVP can only be constructed for well-posed sys- 
tems. While problems in general relativity are typically 
nonlinear we consider here the linear IBVP, as stability 
in the linearized case is a necessary condition for stability 
in the full nonlinear system and these methods may also 
be applied to nonlinear equations. Consider the linear 
IBVP on a domain [0, oo) x Vl 

dtu = A{t,xydiu + B{t,x)u, (1) 
u(0,.x) - f{x), (2) 
w+{t,x) = Sw-{t,x) + g{t,x) ,x dn, (3) 



where u is a vector-valued function, w+ and w- are in- 
coming and outgoing modes, and S is sufficiently small 
(maximally dissipative boundary conditions Q)- The 
system is assumed to be symmetric hyperbolic, i.e., there 
exists a symmetric positive definite matrix H(t,x), the 
symmetrizer, satisfying HA^ — (HA^)^ . The usual proof 
of well-posedness proceeds by defining an "energy" , £ = 
Hu (Px , and by showing that £ can be bounded as 
a function of the initial and boundary data. Analogously, 
a way to construct stable numerical schemes is to design 
them such that energy estimates hold at the discrete level 
This is the procedure we shall follow in this paper. 

Our construction involves four steps. (1) We construct 
discrete derivative operators and a scalar product so that 
a semi-discrete energy estimate holds. The spatial deriva- 
tives are then discretized using these operators, resulting 
in a semi-discrete system of ordinary differential equa- 
tions. (2) We impose boundary conditions in a way that 
does not spoil the previous semi-discrete energy estimate. 
In particular, we account for effects of boundary edges 
and vertices. Steps (1) and (2) guarantee numerical sta- 
bility of the semi-discrete system. (3) We may add artifi- 
cial dissipation and/or arrange the spatial discretization 
to achieve optimal energy bounds. (4) Finally, an ap- 
propriate time integrator is chosen so that fully discrete 
stability holds. 

Step (1): The equations are discretized on a domain 
VL with an inner boundary to accommodate for black 
hole excision. While several grid geometries are pos- 
sible, we concentrate on the simplest case of a cubical 
domain from which an inner, smaller cube has been re- 
moved (both cubes aligned with the grid). We introduce 
the grid points r^fe = {iAx,jAy,kAz) G fl and assume 
that some of these points lie on the boundary dfl. A 
scalar product, S, between any two real vector- valued 
grid functions u and v, is defined as 



N 



AxAyAz ^ aijkufj^Vijk, 

i,j,k=0 



(4) 



where atjk are defined below and in the sum Vijk G A 
semi-discrete energy is then defined hy E — [u,Hu)y.- 

A key ingredient for deriving the continuum energy es- 
timate is integration by parts. Similarly, in order to get 
a semi-discrete energy estimate, one must construct the 
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difference operators, D, and E in such a way that the 
discrete version of integration by parts, caUed summa- 
tion by parts (SBP), holds. In one dimension, SEP is 
expressed by {u,Dv)j: + {Du,v)-s = upfVN — uqVq, and 
this definition can be generalized to higher dimensions 
and more complicated domains. 

A detailed calculation shows that SBP for our chosen 
n holds by defining (similar definitions for y, z directions) 
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Utjk 



Hjk 



at o; = const, faces, 
X — const. OB edges 
and OB vertices 

at a; = const. IB edges 
at IB vertices 
everywhere else 



where IB and OB stand for inner and outer bound- 

'Uijk = {ui+ijk - Uijk) /Ax, 



ary respectively, 

D^^'^Ui-jk = {uijk-Ui^ijk)/ and Dlj^'uijk = (wi+ijA- 

jfe)/(2Ax). The ± signs indicate whether non- 
excised points are to the right or left, respectively, of the 
boundary point. The weights aijk in Q are defined to be 
1 in the interior, 1/2 at the faces of IB and OB, 1/4 at the 
edges of OB, 3/4 at the edges of IB, 1/8 at the vertices 
of OB, 7/8 at the vertices of IB, and zero in the excised 
region. This difference operator is second-order accurate 
at the interior and first-order at boundaries, yielding in 
principle, a scheme with overall second-order accuracy 

Step (2): In contrast to the continuum problem, an 
examination of boundary terms left after SBP in the 
energy estimate indicates that edge and vertex bound- 
ary points make a finite contribution to E at fixed res- 
olution. These contributions show precisely how to im- 
pose boundary conditions at these points. They have 
to be imposed on incoming characteristic modes defined 
by some "effective" normal vectors n. For instance, on 
a uniform grid {Ax = Ay = Az), at an edge one has 
n = {ua + nb)/ ^/2, and at a vertex n — {ria + n})-\- ric) / 
(where Ua, ni, and Hc are unit vectors of the intersecting 
faces). Boundary conditions along these effective direc- 
tions have to be imposed without destroying the semi- 
discrete energy estimate. Following Olsson 55, we do so 
by projecting the right hand side (RHS) of the evolu- 
tion equations at boundary points onto the subspace of 
the grid functions that satisfy the discretizcd boundary 
conditions. The projection is by construction self-adjoint 
with respect to S, which ensures that the solution of the 
projected semi-discrete system will satisfy the boundary 
conditions without compromising numerical stability. 

Step (3): At a fixed resolution, even convergent codes 
may generate significant error growth as time progresses, 
and it is often desirable to minimize this growth. This 
may be achieved by adding artificial dissipation, by rear- 
ranging the semi-discrete equations in a specific manner 
to obtain optimal estimates, or by doing both. We briefiy 
summarize these two procedures. 



Qd^^ + Od""^)"' is sometimes added to the RHS of Eq. QJ, 
which damps high frequency modes but does not affect 
the accuracy of the scheme. Requiring (u, qIj^^u)^ < 
ensures that the semi-discrete energy estimate still holds. 
We have modified the Kreiss-Oliger dissipation operator 
for our black hole excision geometry and S. This modi- 
fied operator has the standard form on the grid interior 
(shown here for the x direction) 



Q(x) 



(where e is a parameter satisfying e > 0.) However, it 

is modified near boundaries to satisfy (u, Q^^^u)^ < 0. 
Let (io,J, fc) label a point on a x = constant boundary 
with neighbors (zq — l,j,k) and (iq + l,j,k). Near this 

(x) 

boundary point becomes 



Q(x) 
d "»o-l 



id "io + l 



(X) T^ix)s 



eAx ,x)2 

Jx)2 



eAx(i^V''^^-2Z?f i^ffV^„+i, 



where the indices jk are suppressed for clarity, and 
is set to zero for points outside the domain. 

A second method for controlling unnecessary growth 
is by rearranging the discretized equations so the op- 
timal estimates for the continuum hold at the semi- 
discrete level. For instance, if at the continuum the sys- 
tem exhibits energy conservation, it is important to pre- 
serve it at the semi-discrete level. For example, assum- 
ing that Eq. has time-independent coefficients and 
that d,{HA') = HB + {HBY, £ is conserved if appro- 
priate boundary conditions are given. The solution of 
dtu — A^{x)DiU + B{x)u may yield growth in E, due to 
the lack of a Leibniz rule at the discrete level. However, 
one can show that rearranging the semi-discrete RHS of 
JTJ as 

gives a non-increasing E, and therefore the discrete so- 
lution will not grow. This idea is closely related to the 
concept of strict stability 5]. 

Step (4): Finally, stability of the fully discrete sys- 
tem follows if it is integrated with a scheme that satisfies 
the (necessary and sufficient) local stability condition or 
the (sufficient) preservation of energy estimate's condi- 
tion (e.g., third or fourth order Runge-Kutta jQ]). 

An example: Wave propagation on a black hole space- 
time g^n presents some of the same challenges as the full 
Einstein equations (e.g. constraint preservation, bound- 
ary conditions, excision of the singularity). Thus it is an 
ideal test bed to demonstrate the techniques described 
in this paper. We write the wave equation in first-order 
form, 



In numerical simulations a dissipative term, [Q 
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V^,, V'T/^^O, V^K=V,1/^, (5) 
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FIG. 1: The self-convergence factor, Q, as a function of time, 

where Q = log2(| |/ai^ — /2Aa; | |/| |/2Aa; — /4Aa: | |), and fAx, f2Ax, 

and f4Ax represent numerical solutions calculated on uniform 
grids with corresponding gridspacing. The inset shows the 
L2 norm of the variable 11 vs. time for the three resolutions. 
The solution is calculated on the domain € [1.5M, 5.5M], 
using 41^, 81^ and 161^ grid points, dissipation parameter 
e = 0, and Courant factor A = 0.8. Since we are outside 
the hole, we chose 6' — 0, giving a symmetric hyperbolic 
formulation. 11 has non-zero initial data of compact support, 
and all other fields are initially set to zero. The solution is 
essentially second order convergent until it decreases several 
orders of magnitude. After that the self-convergence factor 
grows and the evolution is followed roughly until the solution 
for the finest resolution reaches truncation error (about 25 
crossing times). 



with V;^ the covariant derivative. To split these equations 
into evolution equations and constraints, we specify a 
future-directed time-like vector u'^, and contract the first 
and the last equation with it. The evolution equations 



n, 



V'V^ = 0, £uV^ = v^n, (6) 



where £u is the Lie derivative with respect to m''. This 
system is symmetric hyperbolic. Due to the evolution 
equations, the constraints = — = 0, Cf_i„ = 
^ pYv-^ vViL = are Lie-dragged by u^^, i.e., i^„C^ = 0, 
£uCfii/ — 0. 

Although the discussion below applies to any back- 
ground geometry, for definiteness we specialize to the 
Schwarzschild background. This spacetime possesses a 
time-like Killing field, fc'' — (9*)^, and the future di- 
rected unit normal to the t = const, slices is given 
by (fc'' — pf^)/a. We denote with a,/3' and hij the lapse, 
shift and three-metric, respectively. A natural choice for 
the vector field u'^ is n'^, however this requires care spec- 
ifying boundary data compatible with the constraints. 
For a particular coupling between the in- and out-going 
modes, we have found maximally dissipative boundary 
conditions compatible with the constraints. Unfortu- 
nately, they imply reflective boundary conditions in the 
sense that (for homogeneous boundary data and in the 
absence of inner boundaries) the physical energy of the 



FIG. 2: This figure shows the self-convergence factor de- 
fined in Fig. 1, Q. A non-trivial smoothly interpolating is 
used so that the system is everywhere symmetric hyperbolic. 
The domain is x,y,z G [—AM, AM], with an excised region 
x,y,z G [-0.375Af, 0.375M]. The Courant factor and dissi- 
pation are A = 1, e = 0.02, respectively, and runs with 65'', 
129^ and 257^ grid points are used to calculate Q. The initial 
part of the plot shows lower than second order convergence. 
However, experiments in one dimension show a similar effect 
for comparable grid spacings; though the convergence factor 
asymptotically approaches two as resolutions are increased. 



scalar field $ is exactly conserved. Choosing a different 
formulation with = /c^ allows us to impose radia- 
tive boundary conditions, in the sense that the physical 
energy decreases when the wave reaches the boundary. 
In this formulation the constraint variables and C^^ 
propagate tangentially to the boundaries, and thus the 
constraints are automatically satisfied when satisfied ini- 
tially. Moreover, this formulation has the advantage that 
its symmetrizer agrees with the physical energy. Unfor- 
tunately, as the Killing field becomes space-like inside 
the black hole, the symmetrizer is not positive definite in 
this region. 

We combine the advantages of both formulations by 
choosing = {k^ — b^)/a, where is a smooth vector 
field which is tangential to the t = const, slices, vanishes 
in a neighborhood of the outer boundary, and agrees with 
the shift vector as one approaches the event horizon 
from the outside region. This new interpolating formu- 
lation is manifestly symmetric hyperbolic also inside the 
black hole; and at the outer boundary radiative boundary 
conditions can be given. Equations yield 

dtU = b'd,n + —dJam + ^dAVhAm) 



+ -(l3'djb'-b>d,l3')V, + ^d,iVhb')n (7) 
a vh 

dtV^ = d,{aW)+b'd,V, + Vjd,V (8) 

where A* = - b\ H'^ = h'^ - A'A^/a^, and 
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h^^ denotes the inverse three metric. This system is 
symmetric hyperbohc with respect to the energy £ — 
^J^a {IP + H'W.Vj) Vh(Px. The change of £ under 
the flow generated by (|7I8I) can be bounded as 

OB 

where 'w± — fj,±Il + aH'^-' riiVj are the characteristic vari- 
ables which have nonzero speeds, /i± = ±a + /3*ni, and 
rii is the normal to the outer boundary, normalized such 
that K^^riinj — 1. The expression F in Eq. (jSJ vanishes 
if = 0. Then setting w+ to zero guarantees that £ will 
decrease when the wave reaches the boundary. In the pre- 
vious estimate we have assumed that the excised cube is 
sufhciently small (side length smaller than 4\/3M/9) so 
that the inner boundary is purely outflow. 

Equations O-® have been written such that the 
replacement di by Di (only on the dynamical variables) 
automatically leads to a semi-discrete system with energy 
conservation in the case = 0. The semidiscrete solution 
will not grow, even at fixed resolution. 

The figures show the results of the numerical imple- 
mentation of the techniques for the scalar wave example 
we discussed; we have chosen Kerr-Schild coordinates (we 
have checked that the case of Painleve-GuUstrand coor- 



dinates behaves similarly and so does the Maxwell field 
on these backgrounds) . Figure 1 shows evolutions with 
the computational domain outside the event horizon; fig. 
2 shows results of a run with singularity excision. Movies 
can be see in 0. 

In this article we have shown that the use of rigorous 
numerical analysis techniques can be used in practice to 
obtain stable schemes for numerical relativity. This maps 
out a precise road for robust implementations, reducing 
the need for expensive trial and error standard methods. 
The successful application of these techniques to the sim- 
ulation of fields propagating on black hole backgrounds, 
which requires dealing with singularity excision, corners 
and edges, illustrates the potential of the approach. This 
opens the door for applying the same techniques to the 
full equations of general relativity. 
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